function [OD] = getOD(Images) 
%[OD] = getOD(Images)
%Get absorption image from array of images. First layer is atoms, second
%layer is beam, third layer is dark.
%Other common sizes return a stack of OD images.
%Change function so that 1. can accept a stack of [Ny,Nx,3,NImgs] and
%process correctly.
%For different cases, instead of rewriting the code, reformulate so that
%array size is [Ny,Nx,3,NImgs]

%elseif size(Images,3)==5 && size(Images,4)==1
%Images = cat(4,Images(:,:,[2,4,5]),Images(:,:,[3,4,5]));
%elseif size(Images,3)==6 && size(Images,4)==1
%Images = cat(4,Images(:,:,[2,4,6]),Images(:,:,[3,5,6]));
%elseif size(Images,3)==7 && size(Images,4)==1
%Images = cat(4,Images(:,:,[2,5,7]),Images(:,:,[3,6,7]));



if (size(Images,3) == 3)
    atoms = Images(:,:,1);
    beam = Images(:,:,2);
    dark = Images(:,:,3);

    divided = (atoms-dark)./(beam-dark);

    divided(isnan(divided))=1;
    divided(divided<=0)=1;
    divided(isinf(divided))=1;

    OD = -log(divided);
    
elseif (size(Images,3) == 4)
    atoms = Images(:,:,2);
    beam = Images(:,:,3);
    dark = Images(:,:,4);
    
    divided = (atoms-dark)./(beam-dark);

    divided(isnan(divided))=1;
    divided(divided<=0)=1;
    divided(isinf(divided))=1;

    OD = -log(divided);
    
elseif (size(Images,3) == 5)
    atoms1 = Images(:,:,2);
    atoms2 = Images(:,:,3);
    beam1 = Images(:,:,4);
    dark = Images(:,:,5);

    divided1 = (atoms1-dark)./(beam1-dark);
    divided2 = (atoms2-dark)./(beam1-dark);

    divided1(isnan(divided1)|divided1<=0|isinf(divided1))=1;
    divided2(isnan(divided2)|divided2<=0|isinf(divided2))=1;


    OD = cat(3,-log(divided1),-log(divided2));
    
elseif (size(Images,3) == 6)
    atoms1 = Images(:,:,2);
    atoms2 = Images(:,:,3);
    beam1 = Images(:,:,4);
    beam2 = Images(:,:,5);
    dark = Images(:,:,6);

    divided1 = (atoms1-dark)./(beam1-dark);
    divided2 = (atoms2-dark)./(beam2-dark);

    divided1(isnan(divided1)|divided1<=0|isinf(divided1))=1;
    divided2(isnan(divided2)|divided2<=0|isinf(divided2))=1;
    
    OD = cat(3,-log(divided1),-log(divided2));
    %OD = -log(divided1);
    %OD2 = -log(divided2);
elseif (size(Images,3) == 7)
    atoms1 = Images(:,:,2);
    atoms2 = Images(:,:,3);
    beam1 = Images(:,:,5);
    beam2 = Images(:,:,6);
    dark = Images(:,:,7);

    divided1 = (atoms1-dark)./(beam1-dark);
    divided2 = (atoms2-dark)./(beam2-dark);

    divided1(isnan(divided1)|divided1<=0|isinf(divided1))=1;
    divided2(isnan(divided2)|divided2<=0|isinf(divided2))=1;
    
    OD = cat(3,-log(divided1),-log(divided2));
    %OD = -log(divided1);
    %OD2 = -log(divided2);
else
    error('error in getOD: not implemented for %d images',size(Images,3));
end


end

